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INTRODUCTION 


Filtering  of  underwater  acoustic  data  Isfan  operation  which  is 
typically  performed  in  an  analog  manner.  This  operation  is  limited 
in  flexibility  due  to  the  quantity  and  availability  of  hardware  and 
limited  in  accuracy  because  of  component  tolerances,  aging  and  envi- 
ronmental deterioration.  It  is  the  purpose  of  this  memorandum  to 
discuss  the  processing  of  broad  band  acoustic  data  using  digital 
computer  filtering  on  digitized  data.  This  procedure  is  limited  by 
the  need  for  sampling  rates  several  times  the  maximum  frequency  present 
in  the  signal,  but  bypasses  the  above  limitations  by  allowing  the 
simulation  of  a large  number  of  filters.  The  accuracy  is  improved 
and  further  computer  processing  can  be  applied  to  the  results. 
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While  the  computer  programs  are  presently  implemented  to  simulate 
only  Buttervorth  type  low  pass  (LP)  and  band  pass  (BP)  filter,  they 
are  general  enough  to  allow  expansion  to  include  any  filter  which  is 
representable  in  terms  of  a finite  pole-zero  transfer  function  con- 
figuration. The  formulation  presented  here  produces,  under  mild  con- 
straints, digital  filters  with  the  same  impulse  response  as  the 
equivalent  analog  filter. 

THEORY 

Using  z-transform  notation  (ref  a and  b),  a tim6 
x(kT)  , k=0,l,2... 

can  be  expressed  as  a function  of  a complex  variable  z by  the  transforms 

a)  x(z)=  *(Vry*K 

A linear  filtering  operation,  such  as  that  represented  by  the  finite 
difference  equation: 

a*  ••*(**-$  , 

1=0 

can  be  represented  in  z-transform  notation  by  the^equatioq: 

(la)  Xto  , v>vwoi  HUV  ^—~rc 


This  equation  is  analogous  to  the  continuous -time  filtering  operation 
represented  in  Laplace  transform  notation  as: 

(2)  Z(s)=H(s) *X(s) 

Letting  the  z-transform  of  a sampled  time  function  x(t)  be  denoted 
X*(z)  and  letting  the  sampling  operator  on  a Laplace  transform  be 
denoted  by  (x(s))*  (so  that  (X(s))»=X*(z)  ),  we  can  examine  the 
conditions  that  are  sufficient  to  cause  the  sampled  output  of  a filter 
to  be  equal  to  the  output  of  the  equivalent  digital  filter. 

We  will  assume  that  X*(z)  and  H*(z)  are  equal  to  X(z)  and  H(z). 
That  is,  the  z-transforms  used  in  evaluating  the  output  of  the  digital 
filter  were  derived  by  sampling  the  input  to  and  impulse  response  of 
the  analog  filter,  or  by  an  equivalent  method. 
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*•  Z(,)  “d  Z*U)  18  d8t-lTed  11  *-»  ft— 

(3)  Z*(a)=  ^residues  of  7,(8)  at  the  poles  of7,(s). 

I-?,'1 

where  T ^ is  the  sampling  frequency 

<4)  Z*(*)=A-  |Lz  (s. ^ ) , 1. sn,-t-.£y 

Applying  the  sampling  operator  to  eq.(2),  we  have: 

Z*(*)=(  H(s) *X(s)  )* 

Substituting  from  eq. (4) , we  get: 

(4a)  Z*(z)= 

tM,)  **  bm  m — **  ^ 8 

-r(v  ij*  fevSi'f)) . (1. 

Thus,  if  both  X and  H are  bandlimited  to  less  than  half  the  sampling 
frequency,  then  the  sampled  output  of  the  analog  filter  will  be  equal 
to  the  output  of  the  digital  filter  amplified  by  the  sampling  period. 

NUMERICAL  METHOD 

Many  types  of  filters  can  be  represented  in  terms  of  the  pole-zero 
plane*Urati°n  °f  th*ir  transfer  functions  in  the  complex  frequency 

Their  transfer  functions  can  then  be  written: 


H(s)= 


can  then  be  writti 
JT  (s-ri) 

‘fr(a-pj) 

J=1 
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Using  eq. (3)  reduces  (5)  tos 

H*(«)=  K %_ 

jM 

in  the  z-domain.  Letting: 

pjT 

Aj=e  and  Kj= 

we  have  a sum  of  1st  order  terms: 


JS|  (1-Ajz"1) 


The  Kj  and  Aj  will  either  be  real  or  in  complex  conjugate  pairs. 

Each  conjugate  pair  can  be  combined  to  yield  second  order  terms 
of  the  type: 

Cjl  **Cj2  .z-1 

l-BJ1.z“1*Bj2*»“2 

Thus,  H(z)  can  be  represented  as  the  sum  of  first  and  second  order 
terms  of  the  type  seen  in  (la),  representing  a difference  equation 
in  z- transform  notation. 

The  second  order  terms  are  thus  realized  by  the  equations: 

yj(k)  -Bj1.yj(k-l)+Bj2.yj(k-2)=Cji.x(k)-Cj2.x(k-l)  ,k=0,l,2... 

for  as  many  second  order  terms  as  there  are  (j=l,2,...),  while  the  first 
order  terms  are  realized  by  the  equations: 

yjGO  -Bj.yj(k-l)=Cj.x(k)  , k=0,l,2...  , Bj=Aj,  Cj=K 

again  for  as  many  j as  there  are  first  order  terms. 
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The  outputs  yj (k)  are  then  summed  to  give  the  filter  output 
y(k)  (for  each  k).  Note  that  this  is  justified  by  the  linearity  of 
the  z-transformation. 

The  filter  constructed  above  is  called  a recursive  munftrl ggj  filter 
(RNF),  where  the  term  "recursive"  refers  to  the  fact  that  past  outputs 
contribute  algebraically  to  the  present  output. 

PROGRAM  INFORMATION 

The  computations  above  have  been  implemented  in  a series  of 
subroutines  that  are  called  by  the  driver  subroutine  BTRRNF. 

BTRRNF  computes  the  coefficients  of  a specified  filter  and  passes 
them  to  a routine  which  actually  performs  filtering  operations. 


BTRRNF  is  called  by  the  FORTRAN  statement: 


CALL  BTRRNF (NPOLES , ITYPE, CF, BW, T , IFERR) 


VARIABLE  PURPOSE  TYPE  INPUT/OUTPUT 

NPOLES  The  number  of  poles  in  INTEGER  Input 

LP  version,  where  NPOLES*  10 


ITYPE  A switch  which  determines 

whether  a LP  or  BP  filter 
will  be  generated.  If  INTEGER  Input 
ITYPE=0,  then  LP, 

If  ITYPE=1,  then  BP. 

Note  that  shifting  from 
LP  to  BP  doubles  the 
actual  number  of  poles 
(see  below) 

CF  The  center  frequency  in  Hz  REAL  Input 

BW  The  bandwidth  in  Hz  REAL  Input 

T The  sampling  period  in  REAL  Input 

secs 
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IFERR 


An  error  Indicator  I 
If  NPOLES  .GT.  10, 
an  error  exists,  IFERR 
returns  as  0.  If  all 
goes  well,  IFERR  returns 
as  1 


TABLE  1 


INTEGER 


Output 


The  filter  is  specified  in  the  call,  and  if  IFERR  is  returned  as  1,  the 
filter  is  loaded  in  the  filtering  subroutine. 

If  the  filter  type  is  BP  ( ITYPE=l),  a low  pass  filter  with  the 
appropriate  bandwidth  is  set  up  and  then  bands hif ted  to  produce  the 
BP  filter  by  substituting  (CF)  /s'+s'/(cf)  for  s in  the  complex 
frequency  plane.  This  is  the  operation  which  doubles  the  number 
of  poles  as  mentioned  in  table  1.  Note  also  that  the  center 
frequency  is  the  geometric  mean  of  the  3 dB  points,  not  the  algebraic 
mean. 

As  mentioned  above,  a call  to  BTRRNF  computes  the  coefficients 
of  the  appropriate  filter  and  passes  them  to  a subroutine  called 
RNF.  This  is  the  filtering  routine  and  may  be  called  by  using  the 
FORTRAN  statement: 


CALL  RNF(X,I,N) 
VARIABLE  PURPOSE 


The  input  array 
containing  the  data 
to  be  filtered 


The  output  array 
to  contain  the 
filtered  data 


REAL 

array 

dimension 

X(N) 

REAL 

array 

dimension 

Y(N) 


INPUT/OUTPUT 


Input 


Output 
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Number  of  points  to  be 
filtered  and  size  of 
X and  Y arrays 


INTEGER  Input 


Table  II 

For  user's  purposes,  this  subroutine  can  be  treated  as  a black 
box  which  returns  one  output  value  for  each  input  value.  By  its 
nature,  it  "remembers"  previous  inputs,  even  after  it  has  returned 
control  to  the  main  program.  The  user  need  not  filter  all  the  data 
in  a time  series  at  one  time,  if,  for  example,  the  data  are  arranged 
in  several  records  on  a tape.  Consecutive  calls  to  RNF  using  part 

° J;he  3er*es  at  a (in  the  right  order,  of  course)  will  give  the 
results  as  if  all  the  filtering  were  done  in  one  call.  The 
following  calls  give  identical  results: 

CALL  RNF (X,Y, 1000) 


CALL  RNF(X,Y,500) 

CALL  RNF(X(501),  Y(501),  500) 

The  first  is  slightly  more  economical,  but  the  computed  results  are 
uQ6  same • 

. follows  that  if  the  filter  remembers  previous  inputs,  there 

should  be  a way  to  make  it  forget  them,  in  order  to  pass  unrelated 
call^h^0^?  the  aa“®  f^ter.  This  capability  is  provided  in  a 

if FORTRAN, ^ i^B  * fUUy  r6laXed  State*  ma  call» 

CALL  ZERRNF 

>Aieh  zeros  the  memory  elements  of  the  filter. 

“J  aasociated  subroutines  are  located  on  FASTRAND  file 
&1787/LIBRY  and  are  available  for  general  use. 
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SUMMARY 


The  subroutine  BTRRNF  allows  the  simulation  of  Buttervorth  filters 
on  a digital  computer.  The  method  used  is  impulse-invariant  in  the 
sense  that  the  computer  output  is  identical  to  the  sampled  output 
of  the  equivalent  analog  filter.  It  is  expected  that  the  errors 
introduced  by  the  use  of  a finite  number  system  in  the  computer 
are  much  smaller  than  are  the  errors  introduced  in  analog  devices 
by  variabilities  in  component  specification. 


Steven  R.  van  der  Veen 
Physicist 
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APPENDIX  A:  SAMPLE  PROGRAM 

A record  on  FASTRAND  file  SRV/82875  was  known  to  contain  a 
bottom  reflected  signal  embedded  in  noise.  Since  it  was  known 
that  the  noise  level  was  high  outside  the  frequency  band  100-400  Hz, 
it  was  decided  to  pass  the  signal  through  a bandpass  filter  with  that 
pass band.  A program  listing  of  the  sample  program  and  the  plotted 
CALCOMP  output  (Figure  1)  follow,  demonstrating  the  use  of  the 
BTRRNF  subroutine. 


1* 

d* 

o* 

*** 

a* 

o* 

f* 

d* 

9* 

Id* 

AA* 

U* 

Id* 

A4* 

id* 

AO* 

if* 

Ad* 

A9* 


c 

L 

L 

V. 

V. 

C 

C 

c 

c 


c 

c 

c 


“JJwi'*"51®’10®8  


PKUvkaM  to  ulMOuSTKaTE  uloXIAL  FILTERING  ROUTINE  UTKRNF 


THt  PAK^T  Mtv.OHu  ON  FASfRANo  Hut  SRV/d2d7b  WAS 
KNOWN  To  cONlAJN  A I lMt  StKAt'S  wMjCH  WAS  PREDOMINANTLY 
NOASc  OuTbAuL  ThE  KhEUUcNCT  CiANu  iU0-4U0  H2, 

THEKtF OhE  A)  WAS  UScU  M,_Rt  10  DEMONSTRATE  THfc  USE.  OK 
Tht  UANuKwSS  FlcIcR. 


uIMtNSAOu  Xt*04d)  »Y  (2002)  ,XduF<4U9e>) 

INITIALIZE  KL0I1AN6  ROUTINE  ANU  KfAO  XN  RECORD 


CALL.  PL01S<  AooF*hOV6,3) 
CALL  KLOt <0.#U.»200) 
CALL  PLO» <U«»A.b,-d» 
CALL  FAC |ORT , o) 
UIU.OOj 
uo  d ls*,2QU0 
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dO* 

Y(IUY(I-l)+.OOb 

di* 

d 

CONT iNUt 

dd* 

Y(iUOi)so. 

id* 

T (dUUd) =A • 

d4* 

CALL  FTKmN(X» *82675* ) 

dd* 

CALL  FTRAN(d,*SRv' »L) 

d®* 

XF(L)99,99, 

dl* 

XStCTR=X 

dtt* 

CALL  FTRaN(0, *SRV,»XS£ClR»d04d*X, XpEKR) 

d9* 

A 

AF  ( APERR-tl ) 9d»  X » 

dO* 

C 

di* 

C 

SET  UP  M b-PCLt  dp  FlLltR  CENTERED  AT  200  HZ#  dOO  H l 

di* 

c 

WXUE  FOk  A SmMPlANO  HA  It.  OF  bUOO  m2. 

dd* 

c 

d4* 

CALL  BTKKNF<b»X'dOO.,dOU.',OU02'XF£RK) 

dd* 

AF ( AFERR  .HE,  X)  00  TO  9/ 

do* 

c 

d7* 

c 

PLOT  THt  UNFiLTERtU  SXO^AL 

dd* 

c 

d9* 

CALL  AAXd(0«»U.»'S10NAL  MMPLATUOt* ,X6»d. »90. ,-,0b» .05»A0, ) 

40* 

A(20UX)=..Ob 

41* 

A<2U0d)=.Qb 

4d* 

CALL  k.lNt(t’A'2U0O'X'0'U) 

4d* 

CALL  SYMoOLTd.  «d.b.  ,.ib» 'uNFXlTEHEO  UaTA ' . 0 , , x;>) 

44* 

CALL  PLOT <0»»4.»-d) 

43* 

4 

40* 

c 

FILTER  *000  POANIS  OF  ltlE  X ARRmY  A no  PL«T  THE  OUTPUT 

47* 

c 

40* 

CALL  KI^F  (X.A.^UOU) 

4** 

CALL  AXIs(U«»U.»*S10NAE  *MPtiTUDt* ,lb»d. »9U. ,-.0b» .0S*10. ) 

SO* 

CALL  LANE ( 1 * A. 2000 . 1 . 0 » U | 

»a* 

CALL  SYMdOL(d.»d.b.(xb»'FlLTcHEU  U<,T*'.U,»ld| 

Sd* 

CacL  PLOl  Ud.,-4.*-d) 

sd* 

c 

S4* 

STOP 

ad* 

99 

STOP  99 

so* 

9d 

STOP  9o 

! if* 

97 

STOP  97 

so* 

enu 

FILTEREO  OATH 


Ti.-'nre  1:  plot  of  output  fror.  sample 
time  vs.  signal  aiiplit 
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APPENDIX  Bs  PROGRAM  LISTINGS 

All  routines  written  for  thi3  program  were  written  in  FORTRAN  V 
for  the  UNIVAC  1108  with  the  exception  of  the  subroutine  RNF,  which 
was  written  in  1108  Assembly  Language.  Listings  of  all  programs 
follow. 


.»**2£S 


00101 

1* 

SUdR JUT 1 wt  b 1 RRNF ( NPOLES*  J TYPE » CF » BW , T , IFERR ) 

00101 

2* 

C 

ooioi 

3* 

C 

**  SUbROuT  iNE  til  RRNF  EVALUATES  1HE  COEFFICIENTS  Of 

00101 

4* 

C 

*♦  A BuTT£K*ORTH  FILTER.  THE  PARAMETER  I TYPE 

00101 

5* 

C 

♦*  CONTROLS  THE  TYPE  Of  FILTER  ILP  OR  tiPl.  IF  ITYPE 

ooioi 

6* 

C 

**  IS  0'  A LOW  PASS  FILTER  IS  GENERATED.  IF  IT 

00101 

7* 

C 

**  is  i»  a bANU  pass  filter  is  generated*  npoles  is 

00101 

64 

C 

**  THE  NUMtiER  OF  POLES  IN  THE  LP  VERSION.  THIS 

Q0101 

9* 

c 

**  VERSION  IS  bANOSHIFTEO  IF  A BP  FILTER  IS  DESIRED. 

OOIOI 

10* 

c 

**  tiW  is  The  BANDWIDTH  of  the  FILTER#  AND  CF  IS  THE 

ooioi 

11* 

c 

**  CENTER  FREQUENCY  (If  ONE  is  NECESSARY).  T IS  THt 

ooioi 

12* 

c 

**  SAMPLING  INTERVAL.  AlL  ARE  IN  UNITS  OF  HZ  OR 

ooioi 

13* 

c 

*•  SECS*  IFERK  IS  AN  ERROR  INDICATOR  ANO  WILL  BE 

OOIOI 

14* 

c 

**  1 ON  A NORMAL  RETURN.  NPOLES  MUST  BE  ,Lb.  10. 

OOIOI 

15* 

c 

•*  IN  all  cases*  if  iferr  IS  SET  TO  1»  a FILTER 

OOIOI 

16* 

c 

OOIOI 

17* 

c 

**  HAS  BEEN  LOAOEO  IN  THE  RNF  FILTER  SUBROUTINE. 

OOIOI 

10* 

c 

•*  SEE  NUSC/NL  TECH  MEMO  221 1*309-70.  *A  PROGRAM 

ooioi 

19* 

c 

**  TO  SIMULATE  BUTTERWORTH  FILTERS  ON  A DIGITAL 

OOIOI 

20* 

c 

**  computer* 

ooioi 

21* 

c 

00103 

22* 

OOUBlC  PRECISION  POLE(2»5)»POLES<2»10)*ZERO(2*1)* 

00103 

23* 

* ZEROS (2*10) * SCALE  * TC 

■ 

00103 

24* 

c 

< 

00103 

25* 

c 

*•  CONVERT  NPOLES* T to  STORE  IN  PROGRAM 

00103 

26* 

c 

00100 

27* 

TC=T 

00105 

20* 

NsNPuLES 

00105 

00106 

29* 

30* 

c 

Ni=0  NUSC/NL  Tech  Memo 

„ . . „ No.  2211-309-70 

00106 

31* 

c 

**  N MuST  tit  ,Lt,  10.  TEST  It 

00106 

32* 

c 

00107 

33* 

IF (N  *GT.  10)  GO  TO  100 

00107 

34* 

c 

00107 

35* 

c 

*♦  LUCmTE  POLES  OF  N-POLE  LP  FILTER  WITH  BANDWIDTH  BW 

Q0107 

36* 

c 

00111 

37* 

CALL  BTRwKT(N*BW»P0L£*sCALE*M1*M2) 

00111 

30* 

c 

00111 

39* 

c 

•*  IF  1TTPE  .EG.  1*  BANDPASS  FILTER  IS  REGUESTED. 

00111 

4U* 

c 

♦*  generate  it  at  STATEMENT  10 

00111 

41* 

c 

00112 

42* 

IFUfYPE  .NE*  0)  GO  ro  10 

00112 

43* 

c 

00112 

44* 

c 

•*  *ERo  ALL  STORAGE  IN  RnF  TO  PREPARE  IT  TO 

00112 

45* 

c 

**  ACCEPT  NEW  FILTER  CoEFFS 

00112 

46* 

c 

00114 

47* 

CALL  LODRNF 

00114 

40* 

c 

00114 

49* 

c 

**  COMpUTE  COEFFS  OF  FILTER  AND  PASS  THEM  TO  RnF 

00114 

50* 

c 

00115 

51* 

CALL  RECURS ( ZERO  * N) » POLE * Ml » M2 , SC AlE  * TC ) 

00115 

52* 

c 

00115 

53* 

c 

**  ALL  O.K.  RETURN 

00115 

54* 

c 

00116 

55* 

IFERRSI 

00117 

56* 

RtTURN 

IMS  PICE  IS  BEST  QUALITY  PRACTICABLE 
raoi  OQPY  JURAISHSD  TO  DDC 


00117 

00117 

00117 

001*0 

001*0 

001*0 

001*0 

00120 

001*0 

00121 

001*1 

00121 

00121 

00121 

00123 

00123 

00123 

00123 

00124 
00124 
001*4 

00124 

00125 

00126 
00126 


57*  C 

56*  C **  CONvtKT  FROM  LOMPASS  TO  BANDPASS  IN  S-DOMAIN 
59*  C 

60*  10  CAUL  LPTOdp  ( CF  ♦ POLL » Ml » M2  * ZERO,  Nl  # POLES  » M3 » M* • 

61*  ♦2tR0S»N2»SCALfc,IFRR) 

62*  C 

63*  C *♦  TtSl  FOR  ERROR  IN  LPToBP  NUSC/NL  Tech  Memo 

64*  C **  IF  IFRR  IS  0»  AN  ERROR  EXISTS  No.  2211-309-70 

65*  C 

66*  IF (ipRR  .tQ,  0)  60  To  100 

67*  C 

66*  C **  2ER0  ALL  ST0RA6E  IN  RnF  To  PREPARE  IT  TO 

69*  C **  ACCfcPT  NEW  FILTER  CoEFFS 

70*  C 

71*  CAUL  LODRNF 

72*  C 

73*  C ♦*  COMPUTE  COEFFS  OF  FILTER  AND  PASS  THEM  TO  RNF 
74*  C 

75*  CALL  RECURS ( 2ER0S » N2  f POLES » M3 » M4 • SCALE  t TC ) 

76*  C 

77*  C **  ALL  O.K.  RETURN 
76*  C 

79*  IFERRSl 

60*  RETURN 

61*  C 


00126  62* 

00126  63* 

00127  64* 

00130  65* 

•0131  86* 


C **  ERROR  RETURN  ttUH  IFERR  set  TO  0 
C 

100  1FERRS0 

return 

END 
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SUBHUU  u Ht  Ktc  Oh  b ( ZERO  , H , HULL  * Ml , , SCALE  * 1 ) 


* THIS  SUuKuUTiNt  PERFORM^,  THE  1MpULSE-1NVAKIaN1  TRANSFORM  FROM  THE 

* b-HLANfc  TO  IHt  L-PLaNfe.  DOUBLE  PRECISION  lb  REQUIRED  FOR 

* AcCUHACi.  ixo TL  I HA | RECURS  NILl  NOf  HANDLE  MULTIPLE  POLES. 

* THE  ZERuS  OF  THfe  FIlTEK  ARE  TO  bE  N IN  NUMBER.  AND  STORED  IN  ZEKO(I.J) 

* WHERE  ZtHOU.J)  AS  THE  REAL  PAR J yF  THt  U-Tii  ZERO.  AND  ZER0T2.J)  IS  THE 

* IMAGINARY  PmrT  Oh  The  U-TH  ZERO, 

* THE  PGLtS  OF  THE  FILTER  ARE  TO  bE  STORED  IN  POLEU.J)  WHERE  I HAS  THE 

* SAME  MEmNINO  AS  bEFoRE.  lh  J Is  l^SS  THAN  OF  EQUAL  TO  Ml,  THE  POLE  IS  A 

* REAL  POLE  Awo  ITS  IMAGINARY  PART  *lwL  bE  ZEnu.  IF  J IS  GREATER  THEN  Ml 

* THE  POLE  WAlL  bt  ONE  Oh  A CuMPtfcX  PAIR.  ONLy  OWE  OF  A COMPLEX  PAIR  IS 

* INCLUDED  I"  THt  POLt  AKRAY.  THE  OtHER  BEING  PRESENT  BY  IMPLICATION. 

* THERE  ARE  <0  bE  Ml  REAL  POLES  A|«U  M2  CONJUGATE  POLE  PAIRS.  SCALE  IS 

* the  filter  amplification  constant,  and  t is  the  sample  interval,  all 

* NON-INTlGEK  VARIABLES  Are  To  BE  DoUBLE  PRECISION  EXCEPT  FOR  THE 

* ARRAY  CoNTAIjmINu  THt  COmPUTEO  COEFFICIENTS  FOR  THE  RNF  ROUTINE, 

* SEt  NUSc/NL  TECH  MEMO  2211-409-70,  »A  COMPUTER  PROGRAM 

* TO  simulate  bUTTERWORTH  FILTERS  On  A DIGITAL  computer* 


DOUBLE  PkECISION  I.ScALE.TBAR 

double  Precision  Z(2) ,p<n ,temp.*<2) ,pole(2*2U) .zero(2,2o) ,tmp(2) 

REAL  C0EhF(4j 

MT0TsM1+m2 

DO  100  Il=l*MTOT 

Z(U=l.OwO 

P(i)=1.0u0 

P(2)-0,0u0 

4T2)S0.0u0 

•*  TEST  F0«  NO  ^EKOS 
IF (N  ,E«,  0)  GO  TO  HI 

**  COMPUTE  PROuUCT  OF  ( POLE (11) -ZERO ( I ) ) FOR  A^u  I AND  STORE  IN  Z 
DO  HO  Is I.N 

TEMPSZ (1) ♦ (Pole ( 1#  in-ZER0(  1.1)1-  2(2)*lPULE(2,lI>-ZtR0{2,I>  ) 
Z(2)-Z(2)*( PulE (I'll) -ZERO (1,1))*  z t 1 ) *lPULt u, I I ) -ZERO (2, I ) ) 

110  Z(1)sTEMh 


rrr  r r r crc  ccr  r.  rr  rrrrrr  cere 
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noi  OOM  IUMUSHH)  XODDQ 


*•  LO^UIt  TMOewLi  Or  tPOL;-UlJ-POLE<I))  FO*  Acl  I ,Nfc.  II  ANO  STOKE 
♦ ♦ IK  P 

lit  l«.U  I.l.H.Ol 

Ire*  *Lt*.  U>  00  ly  ***# 

ltMP=.Pdj*»F'wLfcd*li)-POetUrin  - Pu)*(P0LlU»11>-P0LEU,IM 
P(c)-F>U,*ePwLtd»ll>-POw.cU*in  ♦ Pd)*(P0LtlS:.ll>-P0LEl2r  J>> 

Pd>-lC»i 
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♦ • COKwOOAit  ►'wet  PAIK.  CwMPUIE  TlKm  FOK  CONJUGATE  OF  POLE(I) 

TtMPiPdl *(pOLtd»lI)-P0etUrl)  > - P(2)*(P0Lt.l2»Il)-FP0LE(2rI) > 
►,(k)=P(2j*ePwLt(I»m-H0eUlrl) ) ♦ P d ) * i POLi.  1 2 » 1 1 ) +POLE  { 2 # I > ) 
FUlslLMp 
130  CONTINUE 
I2u  CONI INUt 

*♦  LOCATE  hOl^  IN  a-PLmNE 

Ad)  = OEaP (►'wet  ( 1 r II ) «T)  «UC0^(H0lE(2»  II  >*i ) 

A(2)=  OEaP(PwlE  d » 11 ) *T)  ♦OSIMPOl*- < 2 1 1 1 ) * 1 ) 

**  COMPUTE  2/p*JsCa*.£.*|,  *MPLIP ICmT ION  CONSTANT 

TE*PsPd)*P<lJ+P{*)*PU> 

T6AH=bCAct/«L,-,P*T 

I HP  Cl)-  |UAK*(p(lJ*2el>^(2)*Z(2j ) 

TMP(2)s  |t3MN*(Pd  W(2Wd)*P<2>  ) 

M TLbl  FOn  KtAL  POLL 

IF ( 11  .Ll.  Ml)  00  TO  ISO 


**  COMPLEX  Poll— btCONu  uk^eh  filtek 
OOLFF  dU2.*lr,PU) 

COEFF  U)=2.*(A(I)*  T WP  ( 1 ) ♦A  (t ) * I He-  ( i> ) > 

COtFF (3) l2, *« ( 1 ) 

COEFF (A)lA(1>«A(1)^A(2)*m(2) 

**  LOAD  COtFFlClENTb 

CALL  OhU«.(tOLFF  (I)  »COtFF  »2>  ryOtFr  t j)  ,eOtF>-  (4) ) 
WU  lo  lOu 

•»  htAL  FOtE--F*KSI  OHOtH  fILUK 

ISO  LOEFF(l)iT«Pd> 

COEFF U)iA(l) 

*•  LCAW  COcFFlcitNTb 

CALC  OKOa (C«tr F (1) , COEFF (2) » 

100  wONllNOL 
Ht  ( OHu 

E.xL  _ 
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QOlOl 

1* 

SUiRoUl  lNc  LPlOtiP<CF,Pl,Ml,M2»Zl»Ni»P2,M3,M4* 

00101 

2* 

» N2  * -aC  * iFiiKh) 

00101 

3* 

C 

00101 

4* 

C 

**  lHlb  SUdKOUIlNE  PERFORMS  THE  BANDSHIFT  NfcCESSrtRT 

00101 

5* 

C 

*♦  TO  CONVcHT  A POlE-ZERO  SET  CENTERED  ON  THE  REAL 

00101 

6* 

C 

♦*  AXlb  TO  A PAIR  OF  SETS  CENTEKEO  AT  +-(CF).  IT 

00101 

7* 

c 

**  uOfcb  THib  at  THE  TRANSFORM  S=SPR1ME/CF+CF/SPR1ME 

• 

00101 

a* 

c 

**  THE  POLES  ANu  ZEROS  10  BE  TRANSFORMED  ARE  IN  Pi 

00101 

9* 

c 

**  AND  21  AND  THE  RESULTS  OF  THE  TRANSFORM  ARE  IN 

00101 

10* 

c 

**  P2  AND  22.  NOTE  THAT  EACH  POLE  PRODUCES  2 

00101 

11* 

c 

**  POLtS  ANU  1 XEROtAT  THE  ORlGiN)*  WITH  THE  OPPOSITE 

00101 

12* 

c 

**  OCCUR ING  FOR  EACH  ZERO. 

00101 

13* 

c 

00101 

14* 

c 

**  SEE  NUSC/NL  TECH  MEMO  2211-309-70,  'A  PROGRAM 

00101 

15* 

c 

♦*  TO  SIMULATE  OUTTcRWOKTH  FILTtRS  ON  A DIGITAL 

00101 

16* 

c 

**  COMPUTER* 

00101 

17* 

c 

00103 

16* 

Double  PRECISION  Pl<2ilO)»P2U*lO>*Zl(2»5)»Z2(2>5), 

00103 

19* 

*A»d»C»U»R»PHi » PI rCFT 

00104 

20* 

PI=3, 14159263356979336800 

00105 

21* 

CFT£CF*Pl*2t0u0  NUSC/NL  Tech  Memo 

No.  2211-309-70 

00106 

22* 

N2=0 

00107 

23* 

M3s0 

00110 

24* 

M4=0 

OOlll 

25* 

IFERn=1 

00111 

26* 

c 

00111 

27* 

c 

**  MM=NUMBEH  OF  POLES 

00111 

26* 

c 

*«  NN=nUMBER  OF  ZEROS 

00111 

29* 

c 

* 

00112 

30* 

MM=M1*2*M2 

00113 

31* 

NN=Nl 

00113 

32* 

c 

00113 

33* 

c 

**  EACH  POLE  PRODUCcS  A ZERO  AT  THE  ORIGIN  ANU  EACH 

00113 

34* 

c 

**  ZERO  A POLE.  THESE  WILL  CANCEL.  COMPUTE  THE 

00113 

35* 

c 

**  DIFFERENCE  IN  NUMBER.  IF  MUcTIPLE  POLES  OCCUR 

00113 

36* 

c 

**  AT  THE  ORIGIN  (IF  MM-NN  ,LT,  -1)»  AN  ERROR  WILL 

00113 

37* 

c 

**  OCCUR  In  RECURS » SO  TEST  MM-nN. 

00113 

36* 

c 

00114 

39* 

ii=mm-nn 

00115 

40* 

IF (11)1*7*5 

00120 

41* 

1 IF (11*1)2*3*3 

00120 

42* 

c 

00120 

43* 

c 

**  TOO  MaNT  ZEKOS 

00120 

44* 

c 

00123 

45* 

2 IFERRSO 

00123 

46* 

c 

00123 

47* 

c 

♦*  Pol  POLES  A1  OKlulN,  EVEN  IF  THERE  ARE  MoRL  Hi  AN 

00123 

46* 

c 

**  ONE 

00123 

49* 

c 

00124 

50* 

3 Du  4 is-l » I 1 *-l 

00127 

51* 

M3=M3M 

00130 

52* 

P2(l#M3)s-l.u-05 

00131 

53* 

4 P2(2.w3)=0,000 

00133 

54* 

GO  To  7 

n&xisw* 


00133 

00133 

00133 

00134 
00*37 

00140 

00141 

00143 

00144 
00147 
00150 
OOlbO 
00150 
00150 

00150 

00151 

00154 

00155 

00156 

00157 
00160 
00161 
00162 

00163 

00164 


ab*  C 

56*  C •*  MURt  HOLfeS  THAN  ZEROS*  PUT  ZEROS  AT  0R16IN 
57*  C 

56*  b 00  6 Islfil 

59*  N2SN2+1 

60*  Z2(1#N2)=0.0J0 

61*  6 Z2(2*N2)=0,0D0 

62*  7 CONTINUE 

63*  IK (Km) 11*11* 

64*  MT0T=M1*M2 

65*  ICNTSM3 

66*  C 

67*  C •*  SOLVE  TRANSFORMATION  EQUATI0N(CF/P2*P2/CF-P1 ll)=0> 
06*  C **  FOR  NEft  POLLS 

69*  C 

70*  00  10  I=1»MT0T 

71*  Aspl(l'I) 

72*  BsPl(2»I) 

73*  Cs2,0U0*A*B 

74*  DsA*A-u*B‘»4#0U0*CFT*CFT 

7b*  PHISOATAN2(C»0)/2,000 

76*  RsOSoHT (OSQRT IC*C+0*D>  > 

77*  CsK*OCOS ( PHI ) /2  # 000 

76*  O=R*O5IN(PHl)/2,0U0 

79*  ICNTslCNT*! 
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00165  60* 

00166  61* 

00166  62* 

00166  83* 

00166  64* 

00167  65* 

00167  66* 

00167  67* 

00167  66* 

00171  69* 

00172  90* 

00172  91* 

00172  92* 

00172  93* 

00172  94* 

00172  9b* 

00173  96* 

00174  97* 

0017b  96* 

00176  99* 

00177  100* 

00200  101* 

00202  102* 

00203  103* 


P2  < 1 * ICNT  > sA/2 . ODO+C 
PZ  (2  *ICN  0=6/2,000*0 

c 

c **  test  fok  piu)  nun-real 

IF (I  .ST,  Ml)  GO  TO  9 
C 

C **  pi (I)  IS  REAt.*  SECOND  ROOT  CONJUGATE  OF  FIRST 
C 

M4SM4M 
60  TO  10 
C 

c **  pi (i)  Complex#  find  second  root#  (roots  due  to 

C **  cONouuATE  OF  Pl(l)  AHE  EOUAL  TO  THE  CONJUGATE  OK 
C **  HOOTS  DUfc  TO  PIU)) 

C 

9  M4SM4+1 
ItNTslCNTFl 
P2( 1* ICNT )SA/2#0D0*C 
P2(2# 1CN1 )=B/2.00u-D 
M4SM4+1 

10  CONTINUE 

11  continue 

IF (NN)21iZ1« 


■“SSI yss***^* 
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1 

1 

% 

• 

1 1 

No.  2211-309-70 

00203 

104* 

c 

00203 

10b* 

C ♦♦  SOLvt  TrUNSt-ORMAflOK  LOUATI0MCF/Z2+<.2/CF-21tI>=0> 

00*03 

106* 

t **  FuR  ZtHOS 

00*03 

107* 

C 

00206 

106* 

Do  20  1=1 »NN 

00211 

1Q9* 

N2SN2+1 

00212 

110* 

A=21U#I> 

00213 

111* 

6=21 (*'I) 

00*14 

112* 

C=*,UUU*A*B 

00*13 

113* 

0=A* A-d*6-4 . UU0*CFT*CF  T 

00*16 

114* 

Prtl=UATAN*(C»L»>/2f  000 

00217 

116* 

RsjSuHT  (ubttRT  (C*C-»D*0) ) 

00*20 

116* 

CsR*UCOS(PHU/2.0u0 

00221 

117* 

0=K*ui»IMHhI>/2.0JO 

00222 

116* 

2*(1,N2)=H/2.0D0*C 

00223 

119* 

Z*C2»w2)=B/2.0D0*U 

00224 

1*0* 

N2=N2+1 

00226 

121* 

22(1'N*)=A/2.0|)0-C 

00226 

122* 

22 ( 2. N2 >=6/2.000-0 

00227 

123* 

*0  CONTINUE 

00231 

124* 

*1  continue 

00232 

126* 

RETURN 

i 

00233 

126* 

END 

• 

/ * 
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OOIUU 

1* 

C 

OOIUU 

d* 

c 

**  IrtMtu  PkuvIUcS  ImMEuIaTE  lIN*S  BETWEEN  RECURS 

oOiou 

3* 

c 

**  hn(j  RNPt  ANu  CAN  BE  oyPASSEU  »ITH  APPROPRIATE 

OOiUU 

H* 

c 

**  lALi.6  IF  0£6iR£0»  tn  PROVIDING  ENTRY  POINTS 

ooxuu 

6* 

C 

**  'OKuliMfd)'  AND  *ORUi:cA#b,C.u)  • ANU  THEN  STORING 

OOiOO 

o* 

C 

**  A t u»  Ct  ANO  0 IN  AHaaYS 

OOIUU 

7* 

C 

ooiuo 

0* 

C 

**  uEE  NUSC/NC  TECH  MEMu  2211-3U9-7U.  »A  PROGRAM 

OOlOu 

** 

c 

• ♦ 10  blROEATt  bUTTtRWoKTH  ULTtRS  ON  A DIGITAL. 

ooiuo 

10* 

c 

**  computer' 

0010b 

11* 

c 

OOliU 

id* 

subroutine  immeu 

00103 

lo* 

entry  okuKAmj) 

ooiuo 

19* 

CauL  ORUERl(AtB) 

00106 

16* 

RETUkN 

00107 

lo* 

entry  oko2<a»o,c»d) 

OOlll 

17* 

Caul  0HUtR2(A#d»C*l>) 

OOU* 

10* 

KeTUkN 

00113 

19* 

End 

1 


r — 1 

. • % 
tol  Ain 

Kl*r  »Ki,F 

\ '■ 

. KN»- 

ib  THt  PROGRAM 

• HA  f ACTUALLY 

IMPLEMENTS  THE  FILTER  GENERATED  By 

. bUUKoUUNt  rtttjKi,. 

AT  CONTAINS  ENTRY  POINTS  FOR  LOADING  COEFFICIENTS* 

. ►OK 

HfctMXXNb  THt  FILTtH.  FOR  ZtROlNG  AlL  COEFFICIENT  STORAGE,  AND  FOR 

. ACiUaLlY  FILTERING 

AN  ARRAY  OF  POINTS. 

• 

. stt 

NUbC/NL  TECH  MfcMO  2 211-309-7U 

. *A  PROGRAM  TO  SIMULATE  mtkp/mt  Tsoh 

. bUTTeR»OKTH  FIlTLRS  ON  A DIGITAL 

COMPUTER* 

No.  2211-309-70 

KLGNmM 

»l0> 

• 

STORE  DATA  ON  EOC  CNTR  o 

XPAST 

RES 

1 

• 

PREVIOUS  VALUE  OF  INPUT 

VY1 

RES 

10 

• 

PREV  V*L  OF  OUTPUT*  2NU  ORO  FILTS 

YY2 

RES 

10 

• 

next  to  last  vae  of  output*  2nd 

t 

order  filts 

YYU 

RES 

10 

• 

PREV  vAE  OF  OUTPUT*  1ST  ORO  FILTS 

\ 

NIOKOL 

RES 

1 

• 

number  of  first  order  filters 

N^OKOfc 

■ RES 

1 

• 

NUMBER  OF  2NU  ORDER  FILTERS 

CC1 

RES 

10 

• 

C -SUB-1  2NU  ORDER  FILTS 

1 

CC2 

RES 

10 

• 

C-SUB-2  2ND  ORDER  FILTS 

Ubl 

RES 

10 

• 

B-SUB-1  2ND  ORO  FjLTS 

1 

bb2 

RES 

10 

• 

b-sub-2  2nd  ord  filts 

! • 

C011 

RES 

10 

• 

C-SUB-l  1ST  ORD  FILTS 

ttbll 

RES 

10 

t 

B-SUB-1  1ST  ORD  FlETS 

SaVHEG 

RES 

4 

• 

TtMP  storage  to  save  registers 

*(1) 

t 

CODE  UNOER  LOC  CNT«  1 

HNF* 

• 

ENTRY  FOR  CALL  RNpU.Y.N) 

' 

i 

sx 

Xl'SAVKtG 

• 

SAVE  REG  Xl 

sx 

X2.SAVREG+1 

• 

SAVE  REG  X2 

sx 

X3.SAVREG+2 

t 

SAVE  REG  X3 

sx 

X4.SAVKEGO 

• 

SAVE  KEG  X** 

EX 

XJ.0.X11 

• 

ADDR ( X ( 1 ) ) INTO  X3 

lxi»au 

XJ.  1 

• 

1 INTO  X3  INCREMENT 

EX 

X4.1.X11 

• 

AUDR(YU))  INTO  X4 

EXI»AO 

X4»  1 

• 

1 INTO  X4  INCREMENT 

E* 

X2»*2»Xll 

t 

N INTO  X2  FOR  EOOplNG 

ANX.XU 

X2.1 

• 

INITIALIZE  X2 

. 1 

00 T LOO 

S* 

AO 

• 

• 

OUTER  LOOP  FOR  N FILTER  OpS 

0 INTO  AO(T-ACCUMdEATOR} 

EX 

X1.N20RUE 

• 

N20R0E  INTO  Xl 

JGO 

X1.E00P2 

t 

TEST  Xl  NOT  0 AND  INIT  FOR  LOOP 

J 

F1RST0 

• 

N20RDE  0.  RUN  1ST  ORD  FILTS 

E00P2 

EA 

A2.CC1.X1 

• 

• 

LOOP  FOR  2ND  ORO  jIlTS 

CCl-Sud-Xl  INTO  A2 

FH 

A2.0.X3 

• 

A2  TIMES  X-SUb-X3  INTO  a2 

SA 

A2.A4 

• 

A2  INI 0 A4 

EA 

A2.CE2.Xl 

• 

CC2-SUU-X1  I'  'TO  A£ 

FH 

A2.XPAST 

• 

A2  TlHES  XPAST  INtO  A2 

FAN 

A4.A2 

• 

A4  MINUS  a 2 INTO  A4 

ea 

A2.YY2.X1 

• 

TV2-SUB-X1  INTO  Ag 

* 

FH 

A2.Ub2.Xi 

• 

A2  TIf.ES  BB2-SUB-X1  INTO  A2 

FAN 

A4.A2 

• 

A4  MINUS  A2  INTO  a** 

EA  Ak#YYl»Xi 

• 

YYI-SoB-XI  INTO  A* 

*! 

SA 

A2.YV2.X1 

• 

A2  INTO  YY2— sUB-Xi 

. 

FH 

A2.0bl.Xi 

• 

A2  TlF.ES  Hbl-SUO-xl  INI o A«? 

FA 

A4.A2 

t 

A4  plus  a2  into  *4 

♦ 

. * 

raoa  cx 

SA 

A4.YVl.Al 

FA 

AU.M4 

oou 

XI .LOOPZ 

FIKSTO 

LX 

Xi.MOKL’E 

JGL 

Xl. LOOPI 

J 

STORE 1 

loopi 

LA 

A2.CCU.Xl 

FM 

A2.0.XO 

SA 

A2.M4 

LA 

A2.BUll.Xl 

FM 

A2.TYll.Xl 

FA 

A4.A2 

SA 

A4.YYll.Xl 

FA 

A0.A4 

J6D 

XI. LOOPI 

STOREY 

LA 

A4,0.*X3 

SA 

A4, XPAST 

SA 

A0,0.*X4 

JGC 

X2.0UTlOO 

LX 

Xl.SAVREG 

LX 

X2. SAVRtG+1 

LX 

XO.SAVRLG+2 

LX 

X4.SAVKEG+3 

J 

4.X11 

LOORNF* 

LX.  Xu 

X12.92 

J 

*♦2 

ZEHRNF* 

LX*  XU 

X12.30 

LOOPZ 

sz 

XPAST .X12 

JGC 

X12. LOoPZ 

J 

l.Xll 

OROER2* 

sx 

Xl.SAVHLG 

LX 

Xl»N20RuE 

LA 

A0,*0.X11 

SA 

A0.CC1.X1 

LA 

A0.*l'Xli 

SA 

A0.CC2.X1 

LA 

A0.*2»X11 

SA 

AO.bbl.Al 

LA 

aO*«3.X11 

SA 

A0.bb2.Xl 

AX.  XU 

xl.l 

SX 

X1.N20ROE 

LX 

Xl.SAVREG 

J 

S»Xll 

uRoEHl* 
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A4  INTO  YY1-SU0-X1 
AO  PLUS  A4  INTO  AO(Y  ACCUMULATOR) 
NEXT  2NU  OKU  FILT » IF  ANY 
START  OF  CODE  FOR  1ST  ORD  FlLTS 
N10RUE  INTO  XI 

TEST  XI  NOT  0 AND  INI!  FOR  LOOP 
NlOROE  0.  STOKE  Y 
LOOP  FOR  1ST  ORU  FlLTS 
CCli-SUU-Xl  INTO  A2 
A2  TIMES  X-S08-X3  INTO  A2 
A2  INTO  A4 
bBli-SUb-Xl  INTO  A2 
A2  TIMES  YV11-SUB-X1  INTO  A2 
A4  PLUS  A2  INTO  A4 
A4  INTO  YY11-SUB-X1 
AO  PLUS  A4  INTO  AO(Y  ACCUMULATOR) 
NEXT  1ST  ORU  FlLTt  IF  ANY 
X-SUb-X3  INTO  A4  THEN  INCREMENT 
X3  FOR  NEXT  X VALUE 
A4  INTO  XPAST 

AO  INTO  Y-SUd-X4  THEN  INCREMENT 
X4  FOR  NEXT  Y 
NEXT  POINT  FOR  FILTERING 
RESTORE  X! 

RESTORE  X2 
RESTORE  X3 
RESTORE  X4 
RETURN 

ENTRY  FOR  CALL  LOORNF 
ZEROS  ALL  STORAGE 
92  INTO  X12 

ENTRY  FOR  CALL  ZERRNF 

ZEROS  ALL  NON-COEFF.  STORAGE 

30  INTO  Xl2 

LOOP  ON  ZEROING  STEP 

0 INTO  XPAST+X12 

ZERO  NEXT  LGC.  IF  ANY 

RETURN 

LNTRY  FOR  CALL  0RDER21A.B.C.0) 
PASSES  2ND  ORD  FlLT  COEFFS 
SAVE  HEto  XI 
i«2orde  into  xi 
A INTO  AO 

AO  INTO  Ctl-SUB-Xl 

b INTO  AO 

AO  INTO  CCZ-SU8-X1 

C INTO  AO 

AO  INTO  BB1-SUB-X1 

U INTO  AO 

AO  INTO  HB2-SUB-X1 

INCREMENT  XI 

STORE  NEM  NZOROL 

RESTORE  Xt 

RETURN 

LNTRY  FOR  CALl  0RuF.h1<A,B) 

PASStS  1ST  OK0  FILT  COf.FFS 
SAVE  REG  XI 


Xl.SAVHLG 


LX 

Xl.NluRjt 

LA 

AO.AU'aU 

SA 

AO.CCll'Xl 

LA 

AO,*I»Ail 

SA 

AO.udll'Xl 

AX  »Xu 

Xi.I 

sx 

Xl.NlOKUt* 

LX 

Xl.SAVwtG 

J 

END 

J'Xil 
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NioRUe  into  xi 
a into  ao 

AO  INTO  Crli-SUb-Xl 
0 INTO  AO 

AO  INTO  BRll-SUtt-Xl 
INCREMENT  XI 
STORE  NE*  NiORDE 
RESTORE  XI 
RETURN 


